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Abstract 

We formulate quantum optics to include frequency dependence in the modeling of optical networks. 
Entangled light pulses available for quantum cryptography are entangled not only in polarization but 
also, whether one wants it or not, in frequency. We model effects of the frequency spectrum of faint 
polarization-entangled light pulses on detection statistics. For instance, we show how polarization en- 
tanglement combines with frequency entanglement in the variation of detection statistics with pulse 
energy. 

Attention is paid not only to single-photon light states but also to multi-photon states. These are 
needed (1) to analyze the dependence of statistics on energy and (2) to help in calibrating fiber couplers, 
lasers and other devices, even when their desired use is for the generation of single-photon light. 

PACS numbers: 03.65.-w, 03.65.Nk, 03.65.Ta, 84.30.Sk 
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PART II 



10. MODELING POLARIZATION-ENTANGLED QKD 

The language offered in the preceding sections supports a wide variety of models of the 
system shown in Fig.|5]|l5^], as well as of many other systems. Here we offer a first-cut model 
of a fiber-optic network that employs polarization-entangled light for quantum key distribution 
(QKD). At the present stage of development of quantum key distribution, the purpose of the 
modeling can hardly be to replace experiments: we lack convincing reasons, whether theoreti- 
cal or experimental, on which to ground the guesswork necessary to generate numbers. Rather, 
we show how a model drawing on some questionable guesses can stimulate experiments. 

Specifically, the model offered shows that if one assumes a Poisson distribution and one 
assumes a light state invariant under identical SU(2) transformations of the light to both Alice 
and Bob, then such and such relations hold between energy, Bob's error rate, and Evangeline's 
entropy. These relations are the "conclusion part" of a statement that includes also an "if 
part." As remarked by Dave Pearson, any interesting conclusion makes the "if part" worth 
exploring, for instance by means of experiments. What evidence can we find in the lab for or 
against the assumption of SU(2) invariance? For or against a Poisson distribution of photon 
number? Expecting to later challenge some of our own assumptions, we try to make our 
modeling modular, so that the assumptions can be changed, both to make room in the future 
for improvement, and to lay the ground for studies of sensitivity to these assumptions. 

Now to business. Consider a simplified experiment to explore the relation between the 
energy of a polarization-entangled light pulse and various detection probabilities relevant to: 
(1) quantum bit error rate (QBER), (2) sifted bit rate, and (3) an eavesdropper's entropy. We 
start by being interested in the QKD topology shown in Fig. 13 in which Alice and Bob each 
operate passively to detect in two bases, with two detectors per basis. 

For this first model we assume: 

1. An ensemble of trials, one entangled pulse pair per trial. 
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FIG. 5: Polarization-entangled QKD system subjected to eavesdropping attack. 

2. Pulses timed and detectors managed so that trapping, dead time, and other memory ef- 
fects of detection are negligible. 

3. Light state leaving transmitter, propagating in modes oi, 02 to Alice and in modes bi, 62 
toward Bob, invariant under application of any given SU(2) transformation to both the 
a-modes and the 6-modes, as discussed below. 

4. A probability of 1/2 for Alice and Bob having matching bases. 

5. Alice and Bob use "on-off" detectors as described in Sec. [HI the detector for mode aj has 
a dark-count probability Pdark(aj) and an efficiency 77det(ai)- 

6. A fraction r/trans of the energy transmitted to Bob survives attenuation, as described by 
frequency-independent coupling to an undetected (undesired) mode. 

7. Poisson distribution of photon number in the energy transmitted to Bob. 

8. Sifting rule: Alice and Bob discard a bit except when (a) they each get one and only one 
detection, and (b) their bases match. 
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9. Eavesdropping attack: Evangeline sneaks a non-polarizing beam splitter into Bob's 
fiber to siphon off her choice of a fraction of the light into modes 63 and 64 that propagate 
in her fiber. This runs through a long, lossless, delay line to a rapidly variable polarization 
rotator, followed by a pair of perfect photon-number detectors. The delay line allows 
Evangeline to postpone detection until she has learned what basis Bob has used for his 
detection; she then rotates or not, as necessary to choose the basis that matches Bob's. 

Much of the analysis is independent of assumptions 5 and 6; we indicate later where these 
assumptions enter. 

To analyze the assumed eavesdropping attack, in which our attacker Evangeline knows 
Bob's basis, we do not need the whole setup of Fig. |5l we can simplify by leaving out the 
rotated bases, and the two beam splitters that support them, along with Evangeline's variable 
polarization rotator. This results in Fig.|6l One can think of the modes with subscripts 1 and 
3 as 'vertically polarized' and those with subscripts 2 and 4 as 'horizontally polarized.' We 
model the variable coupler by which Evangeline taps off energy by an SU(2) transformation. 
For j = 1, 2 we have 



bj{uj) 




U —V* 




h 








V u* 









(10.1) 



where the are vacuum modes assumed unexcited, and + = 1. Inverting this equation, 
we find 

hj{uj) = u*bj{uj) + v*bj+2{^). (10.2) 
A. Outcomes and probabilities 

By an elementary outcome we mean a possible joint response of all the detectors involved. 
We view an elementary outcome as constituted from components, one component for each 
detector. In the context of the model presented here, an elementary outcome consists of a 
bit string for the modes for which detection is binary (such as APD detectors as modeled in 
Sec. El), and a non-negative integer for each of Evangeline's two modes (that we imagine as 
having photon-counting detectors). Each possible joint response of the binary detectors can be 
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FIG. 6: Simplified network. 

expressed by partitioning the set of modes subject to binary detection into a set Jq for which 
the response is 'no-detect' and a set Ji for which the response is 'detect'; correspondingly 
any elementary outcome has the form (Jq, Ji, A;, m) where k = N{bs) and m = N{b4). The 
corresponding detection operator factors; given any normalized state vector \ip). 



Pr(Jo, Ji, A;, m) 



|M(Jo,Ji,fc,m)|V'), 



(10.3) 



with 

M(Jo, Ji, k, m) = PMPM I n Mo{x) ](l[M,{x)Y (10.4) 

\x£Jo / \xeJi / 

where Pkib-^) and Pm{bi) are projections and x ranges over modes in the lists Jq and Ji. 

When we want to ignore Evangeline's detections, we have a non-elementary outcome 
(Jo, Ji) with the corresponding probability 



Pr(Jo,Ji) = Pr(Jo,Ji,A;,m). 



(10.5) 



fc,m=0 



In expressions for outcomes that enter calculations we will often write, in place of (Jq, Ji), a 
list of all the modes in these sets in the order (ai 02 hi 62) with a bar placed over the (undetected) 
modes that belong to Jq. Thus an outcome specified by Jq = {02, foi} and Ji = {oi, 62} will 
also be written iai 02 &i &2)- 
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Only four such outcomes survive the assumed sifting rule that requires exactly one detection 
by Alice and exactly one detection by Bob. Because of the assumed light state. Bob takes 
'detect' on bi to be a 1-bit for a quantum key, while Alice takes 'detect' not on ai but on 02 for 
a 1-bit, etc. With this rule the four outcomes that survive sifting are related to 'correct bit' and 
'error' as follows: (oi 02 61 62) and (oi 02 &i 62) are correct from the standpoint of QKD, while 
{cLi a2 bi 62) and (ai 02 61 62) are errors. 

The probability of a bit surviving sifting for cases in which the bases match is, in this 
simplified model, 

Pr (trial produces a sifted bit) 

= Pr{aia2bib2) + Pr{aia2bib2) +Pr{aia2bib2) + Pr{aia2bib2). (10.6) 

The probability of an error in a sifted bit is 

o /. I -ft ^^ Pr(aia2 6i62) + Pr(aia2 6i62) 

Pr (bit error sifted) = — —. ^ ^. (10.7) 

Pr(tnal produces a sifted bit) 

For trials that result in error-free bits, we want to know the degree to which Evangeline's 
outcomes N{bs) and N{b4) leave her ignorant concerning Bob's outcomes. This ignorance 
of Evangeline with respect to error-free bits is measured by Renyi entropy of order R. This 
i?-entropy depends on the conditional probability that Bob received a 1, given Evangeline's 
detector response {k, m). This we denote 

Ev(/c,m) =^ Vr{^ia2bib2,k,m\{aia2bib2]k,m} ov {aia2\b2]k,m}) 



Pr(ai 02 61 b2]k,m) 



(10.8) 



Pr(ai a2 bi 62; k, m) + Pr(ai 02 &i 62; k, m) 
Evangeline's i?-entropy given N{b'i) = n and N{bA) = mis 

Bntnik, m) = -^—[{Ew{k, m))^ + (1 - Ev(A;, m))^]. (10.9) 
1 — R 

Evangeline's average /?-entropy on error-free bits is then 

Srm=o[P^(^i "^2 h 62; k, m) + Pr(ai 02 bi 63; k, m)]EntR(A;, m) 

AvEntij = — ^- — . (10.10) 

Pr(ai 02 bi 62) + Pr(ai 02 61 62) 
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Altogether there are six types of outcomes: those with and without the distinction "/c, m" for 
Evangeline's detectors, and for each of these the general case of an arbitrary energy distribution 
and two special cases of the Poisson energy distribution and that of a single photon number n. 
Of course all these probabilities (and the i?-entropy) depend on both the light state l^p) and the 
parameters of the APD model of Alice's and Bob's detectors. 

B. Light state 

We formulate a family of states for entangled light, as is discussed in more detail in Ap- 
pendixO The calculations are complicated; here we carry them out for two limiting cases that 
are relatively simpler. 

To begin rather generally, we are concerned with an otherwise arbitrary normalized state 
vector 

oo 

\^P)=Y,Cn\^n), (10.11) 
n=0 

where the state \ijjn) signifies n photons transmitted to Bob and 

oo 

j]ic„r = i. (10.12) 

n=0 

The index n for 'photon number' to Bob (which here will be the same as that for Alice) takes 
specific meaning when we assume a polarization-entangled light pulse invariant under match- 
ing SU(2) transforms of both a-modes and 6-modes, for which 

\^n) = m9on)[g^:{aib2 - a2biyr\0), (10.13) 

with 

/oo _ _ 

duJiduj2g(;{uJi,uj2)[a\{uJi)bl{uj2) - al{uji)b\{uj2)], (10.14) 
-oo 

where we assume a family of functions g(^(uj,Lj) of the following form. For any real- valued 
functions (p{uj) and 0(cj) and positive real parameters a and a, let 

^C(-,-) = -Le^^(-)e^^(-)Ffc;^^,^) , (10.15) 
Vaa V or a J 
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where we define 



(10.16) 



regardless of the value of C, 

/oo 
dxdi/|F(C;x,i/)|2 = l. (10.17) 
-oo 

Thus for any choice of center frequencies ujq and Cjq, bandwidth parameters cr and cr, and phase 
functions 0(u;) and 4>{Cj), we get a family of (7^'s. 

In Eq. (110.131) . ^{gc^-, n) is a normalization constant that makes |V^„) have unit norm, so that 
it is defined by 

[Ar(^7c,ri)]' = (0|[(7c*:(ai&2-a26i)]"[(7c:(ai62-a26i)t]"|0)-\ (10.18) 
Writing {ajhk) for g'^ : ajfe^, we have from Eq. dlO.lSI) 

W{gon)f = (0|[(ai62)-(a2&i)]"[(ai62)-Mi)]^"|0)~' 

1 (0|(«i&2)'(a2&i)"-'=(ai62)^'(a26i)«"-'=)|0) 



n 



where we define 



dcf 



(0|(^7c*:afer(^7c:at6t)«|0), 



(10.19) 



(10.20) 



and the last equality in Eq. (110.191) comes from tensor-product factoring. Note that all that 
matters about a{uj) and h{uj) in this definition is that they are mutually orthogonal; any other 
pair would give the same value. 

Re-expressing \ipn) in terms of detector modes per Eq. (110.21) . one has for the light state 



IV'n) = M{gi^,n)[g^:{ua\h\ + mj^l - ualh\ - va\hl)f\Q) 



(10.21) 
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When we write {Ojbl) as shorthand for : aj^^,, this becomes 

jlkW.ml 



j,k,l,m=0 
j+k+t+m=n 

(10.22) 



C. Energy profile 

We suppose that the light is generated by equipment close to Alice, so that the energy 
exposed to eavesdropping is in the 6-modes rather than in the a-modes. For a state of the form 
defined by Eqs. (110.1 IL (110. 13I) . we want to express the expectation energy for modes hi and 
1)2, denoted 

Energy(6i,62) = {i^\H{hiMm , (10.23) 
where for present purposes we approximate the hamiltonian operator 

duo \uo\[h{{uj)hi{uo) + hl{uj)h2{uj)\ (10.24) 

-oo 

for narrow-band signals by 

/oo ^ _ _ 

duo [h{{uj)hi{uj) + hl{uj)h2{uj)l (10.25) 
-oo 

where uoq is the carrier angular frequency, similar to that of Eq. (12.281) . The commutation 
relation of Lemma (IB 141) of Appendix |B] leads then to 

oo 

Energy(6i,62) = ^oJ^^IC^I^ (10.26) 
Denote the 'mean photon number' by 



n=l 



/i = -^Energy(6i,62), (10.27) 

so that, from Eq. (110.261) we have 

oo 

^ = ^n\Cn\'^. (10.28) 

n=l 
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To produce a dependence of probabilities on we have to choose an energy profile, which 
means choosing the C^. As a first cut, we will show consequences of assuming a Poisson 
distribution 



D. Calculation of probabilities 

Lacking strong theoretical or experimental evidence to guide the choice of energy distribu- 
tion of light for QKD, we arrange the modeling so that this distribution can be entered as a 
parameter. To this end we provide for modeling the contribution of individual values of photon 
number n. As remarked earlier, we want probabilities Pr(Jo, Ji, k, m) and Pr(Jo, Ji), i.e., 
with and without the "k, m" distinction; further we want each of these for the general case of 
an energy distribution C = {|C„p} and for the two special cases of (1) a Poisson distribution 
and (2) an n-photon state. Altogether, this makes 2x3 = 6 types. Each of these six types will 
be expressed by a corresponding function T; we will soon see Tc^km', %i.,km', ^.fem^ ^c, 7^, and 
Tji. These T (for "total") functions will be calculated as sums of corresponding functions 
that are decorated with the same subscripts. 

To start with, for purposes of calculating probabilities we break the state {i/jn) down fur- 
ther in terms of the response of Evangeline's photon-number detectors, assumed expressed by 
projection operators for modes 63 and 64: 



Assume: |C, 



n 



2 




(10.29) 



n 




(10.30) 



k,m=0 
k+m<n 



where 




dcf 



Pk{h)PUb4)\i^n) 




j=0 



x(-l) 



{a{bly{a{binalb{) 



n—k—m—j 



{albino). 



(10.31) 
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and \il)n) but not \il)n,km) have unit norm 

(^n,fcm|V'n,fcm) = (V'nl^n) = 1- (10.32) 



k,m=0 

As shorthand for products of detection operators we write 



Mo (Jo) = X{M, 



Mi(Ji) = W Mi(x). (10.33) 

The detection operators for Alice and Bob assumed here "respect" the numbers n, k, and m 
in the sense that, for any generic M of these operators, 

oo 

n=0 
n 

(^„|M|^„) = ^ {i'n,km\M\i,n,km). (10.34) 

fc,m=0 
fc+m<n 

Thus from Eq. ( 110.31) we get, putting all this together, 

oo 

Pr(Jo,Ji,fc,m) = ^ |C„|2(^„,fc^|Mo(Jo)Mi(Ji)|V'„,U, (10-35) 

n=k+m 

where we adopt the convention that \il)n,km) = if A; + m > n. For probabilities that are 
indifferent to Evangeline's outcome components, we have 

oo n 

Pr(Jo,Ji) = 5^|C„p (^n,fe,n|Mo(Jo)Mi(Ji)|^„,fc„). (10.36) 

n=0 k,m=0 
k+m<n 

This calculation is centered on 

|Mo(Jo)Mi(Ji)|^„,fc„). (10.37) 

From Eq. (18.231) we have 

(V'n,fcm|Mo(Jo)Mi(Ji)|V'n,fcm) = ^,fcm(Jo, Jl), (10.38) 

where we define 

r„,fcm(Jo, Ji) = -^".fc-lJollX), (10.39) 

XcJi 
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with Tn,km defined by 

^„,,^(Jo||X) (-l)#(-'°ll^)(^„,fc„,|Mo(Jo||X)|^„,fc„,)- (10.40) 

(Note that the sum is over all subsets of Ji, including both Ji itself and the empty set 0.) 

In our numerical programs it proved convenient to code the arguments jF(Jo||X) (of any 
type of jF-function) by a 4-bit vector ordered by all the modes ai, 02, 61, &2, with a 1 if the 
mode belongs to Jo||X and zero otherwise. Thus jF(ai, 027^1) would be coded as JF(IIIO). 
For example, for one of the probabilities that enter Eq. (110.81) . we have 

Tn,km{ai a2 &162) 

= ^ n,,fcm(Ol, ^2) + ^ ra,fcm(Ol, &2) '^2) + ^ ra,fcm(Ol, &2; &l) + ^ n,km{(^i, &2; (12, ^2) 

= Tn,km{l^m + ^n,fc™(1101) + J-„,fc„(1011) + (10.41) 

This is convenient because of a trick of using a second coding scheme for coding the argument 
of the T-function: The T-function argument is coded ('negatively' so to speak) by assigning a 
1 if the mode appears with a bar over it and a otherwise. With these two coding schemes, the 
code for the T-function argument becomes the code for the first jF-function argument, and the 
rest of the .F-function arguments are obtained by "filling in zeros" in all possible ways. 

In this way one evaluates Eq. (110.371) using only Mq operators. Drawing on the prescription 
of Proposition (I8.16I) . to calculate J-'n,km we define 

|0n,fcm(ai,a2, A,/52)), (10.42) 

the vector obtained from the expression in Eq. (110.311) for \ipn,km) by replacing, for j = 1, 2, 
a](co') by ay^a](t<;) and fe](u;) by /5j^^6](co'). This substitution yields 

(0n,fcm(ai, "2, A, /52) |0n,fcm(ai, "2, A, f32)) 



2\k+in 



E "" / j+k n—k—jnn — k—m—jnj 

— ^— -— ^ a, ^dix. 

(10.43) 

where we define 

X = {0\{aihyia,h)\a2hr-'^-'^-^{a2hr\ialbl^^ 
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{Q\{aihy{a,h)''{a\bly{a\bl)''\0) (0|(a26i)"-'-'"-^(a263)'"(4&I)"-'-'"-^'(44)"|0), 

(10.44) 



with the factorization coming from a tensor product. From Proposition (IB 181) we have 

maM{a,b,)\a\bly{a\bl)''\0) = _Zl^(0|(a6)^+^(at6t)i+fc|0) = + fc) 



(j + ky. 



where, as in Appendix O we define 



(10.45) 



(10.46) 



As a result, we have 



X = jl^K"- -k-m- j)!m! Sg (j + k)Eg Jn - k - j) 



(10.47) 



whence it follows that 



(0n,fcm(ai, "2, /32) |0n,fcm(ai, ^2, A, /52)) 



|2(n-fc-m)^]^ _ |^|2^ 



2\A:+m 



k\m\ 



j\{n — k — m — jy. 



(10.48) 



where we define 



Gn,kmiw,X,y,z) 
def 



Wgony'nl' Eg^jj + k)Eg^{n - k ~ j) ^kyU^k-m-j 



k\m\ 



j=0 



j\(n — k ~ m — jy. 



with 



(10.49) 



w = aiP2\u\'^, a; = — |up), y = a2Pi\u\'^, 2; = 0:2(1 — I^P). (10.50) 



Thus for L C {ai, 02, 61, 62}, the recipe of Sec.[Hlyields for the J^n,km of Eq. (110.401) 



:Fn,kM = [nil -Pdark(a;)] ) g^n, 

VxeL / 



(10.51) 
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with Qn,km{w., X, y, z) evaluated, using Eq. (I10.50I) . according to: 



1 - ?7det(aj), if ctj e L, 
1, otherwise, 



1 - iff^eL, ^^^ ^^^ 

1, Otherwise. 

The remaining probabilities that we need to evaluate are less fine-grained; they are 

(V'„|M(ai,62;a2,&l)|^n) = ^ (V'n,fcm|M(ai, 62; 02, &l)|^n,fcm), (10.53) 

k,m=0 
k+m<n 

along with (?/^„|M(a2, 61; ai, &2)|'^n) and the two error outcomes {^pn\^{a2,b2', o-i, bi)\'4^n) and 
('^„|M(ai, 61; 02, &2)|V'n)- These are calculated by replacing Tn^km by 

r„(Jo,Ji)= 5^ r„,fc^(Jo,Ji), (10.54) 

/c,m=0 
A;+m<n 

the calculation of which is streamlined by noticing, in analogy to Eq. ( 110.391 ). 

r„(Jo, Ji) = J2 (-l)*^''^^n(Jo||X), (10.55) 



XcJi 



with 



J^n{'L)= ^n,kM, (10.56) 



k,m=0 
k+m<n 



evaluated using 



gn{w,x,y,z) = ^ Qn,km{w,x,y,z). (10.57) 

fc,m=0 
k+m<n 

Letting r = j + k and re-ordering the sum and using Eq. ( 110.491) . one finds 

gn{w,x,y,z) 



w^x"^ ^ ^— ^ y^z^ ^ 



j\(r-j)\ ^ £Un-r-i)\ 



" ^ (r)'Er,An — r) 



lATig^, n)\'n\' ^ \ w + xJiy + zy-^ (10.58) 

^-^ r\[n — r)\ 



r=0 
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Turning to the probabilities that require summing over a distribution of energies, we use Eq. 
(110.3811 to make Eq. (110.351) explicit: 



oo 

dcf 



FT{Jo,Jl,k,m) =Tc,km = ^ \Cn\'^%i,km, (10.59) 

n=k+m 

where we recognize that Tn,km and Tn^km are zero if < k + m. Similarly Eq. (110.361) becomes 



oo 

def 



Pr(Jo,Ji) = Tc(Jo,Ji) = Yl |Cn|'7;(Jo,Ji). (10.60) 

n=k+m 

Given an energy profile C =^ one evaluates Eq. (110.351 ) efficiently by defining 

oo 

Qc,kmiw,x,y,z) = ^ \Cn\'^gn,kmiw,x,y,z), (10.61) 

n=k+m 

and observing that, analogous to Eq. (110.511) . 



^ C,fcm(L) 



.1)#(L) -pdark(x)] J 6^C,fc,n(w^,a;,l/,2:) bval , (10.62) 



evaluated by the prescription of Eq. (110.521) . 

For evaluating Pr(Jo, Ji) (for use when one is indifferent to Evangeline's outcome compo- 
nents), we introduce the analogous 

J-c(L) = (-1)#(^) (ntl -Pdark(a:)] J gc{w,x,y,z) Invai , (10.63) 

evaluated by the prescription of Eq. (110.521) . but with 

oo 

Qciu],x,y,z)= ^ \Cn\^Gniw,x,y,z). (10.64) 

n=k+m 

E. Case I: No frequency entanglement 

The absence of frequency entanglement is exemplified by gi{uj, cD) = f{uj)h(uj), normalized 
so that J J du duj \g{uj,uj)\'^ = j duj \ f{oj)\'^ = J dCu \h{Cu)\'^ = 1. Denoting n) evaluated 

at C = by Mi{n), we have from the rules for Case I at the end of Appendix O applied to Eq. 
(fTaT9b 

mn)\^=(nl±(Aklin-k)!\ = ^^^^^ (10.65) 



fe=0 
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From Lemma (IC72I) of Appendix O we have 

Si(n) 1 (0| (a ■ ■ 6)"(at ■ ■ 6t)"|o) = n\ . (10.66) 
With these specializations, Eq. (I10.49I) becomes 

= I ^fe^n. "y'" ij + k)\{n-k- j)\ ^,yn-k-m-,_ ^^^^^^ 

(n + l)A;!m! j\{n — k — m ~ j)\ 

« — I 

For reference, we note that this involves a hypergeometric function [.16.] 

X V"''"'"^'" 2i^i ( A; + 1, + m - n; - n; - ) . (10.68) 



{n + l)m\{n — k — m)\ \ ' ' ' 1/ 

For use when one is indifferent to Evangeline's outcome components, one finds {w, x, y, z 
analytically from Eq. ( 110.581 ) as 

'(I)^„ ^ „ def (I) 



A:,m=0 
A;+m<n 

{y + ^)"+i - (w + x)"+i 



(10.69) 



[n + l){y + z — w — x) 
In numerical work, we encounter the limit as y + z — w — x ^ 0, in which case this becomes 

lim {w, X, y, z) = {y + z^ . (10.70) 

y+z—w—x^^ 

For the special case of the Poisson energy distribution, we sum Eq. (110.691 ) to obtain 

oo „ 

gf{w,x,y,z) =^ e-^^ — Q''^^[w,x,y,z) 



n=0 



- [e'^{j/+-) - e^^'^+^^j . (10.71) 



y + z — w — X jj, 

In numerical work, we again encounter the limit d&y + z — w — x 0, in which case this 



becomes 



lim g^^\w,x,y,z)=e-^''^^-y-'\ (10.72) 



57 



F. Case II: Limit of extreme frequency entanglement as C ^ ±00 



It makes no sense to ask for the limit of gc_; however, we explore how the probabilities 
behave in the limit of large values of \C\. Denoting Af{g(, n) in the limit as C ^ ±00 by 
A/ii(n), we have from the rules for Case II at the end of Appendix ICl applied to Eq. (110.191 ) 

2-n 



IMi(n) 



(10.73) 



which (in its power of n\) differs from the preceding case. From Lemma (IC79I) of Appendix ICl 
we have 



Sn(n) '=2? lim 1 (0|(a • C ■ hj [a^ ■ ^ ■ = 1- 



C-^±oo n\ 

With these specializations, Eq. ( 110.491 ) becomes, for Case II 



(10.74) 



^n,L(^'3;,y,z) 



2-"n! 



n—k—m 



k\m\ j\{n — k — m — i] 



2-"n! 



n—k—m 



k\m\{n — k — m)\ 

Putting this together with the Poisson distribution for energy yields 

00 



(10.75) 



n=k+m 



E 

n=k+m 



/i\ " X z'^{w + y) 



n—k—m 



k\m\{n — k — m)\ 



For the sum of these over k, m, Eqs. (110.741) and (110.581) imply 



k ^m 



k\m\ 



2-"(w + x + i/ + 2)". 



(10.76) 



(10.77) 



[w,x,y,z) 

Summing over all n weighted by |C„p yields, for use in calculating probabilities for Alice and 
Bob regardless of Evangeline's outcome components, 

def 



Gi^^\w,x,y,z) 



' J2 ^2'''{w + x + y + zr 

n=k+m 

exp |^~^(2 — w — X — y — z) . 



(10.78) 
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G. Example numbers 

Figure III shows the probability of error vs. fj, for several values of C- This shows explicitly 
how changing frequency entanglement C, changes the dependence of the probability of error on 
ji. Other cases can be generated from the MATLAB programs listed in Appendix 10 

APPENDIX A: BACKGROUND 
1. Commutation relation 

The commutation relation was chosen by analogy with that for the harmonic oscillator. One 
can ask if the (5-function should be multiplied by a factor that depends on propagation constant. 
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We answer "no" for the following reason. We want the energy of a 1 -photon light state having 
a narrow frequency spectrum centered at uq to be close to h\tuo\. Taking such a state as an 
in-state to a fiber-vacuum interface results in an out-state of the same frequency but different 
wavelength. In order for energy to be conserved, we need the energy calculated for a 1 -photon 
state to be independent of variations in wavelength at a given frequency. That rules out any 
factor in the commutation relation that depends on the in-fiber propagation constant. 

2. Units 

a{Lu) in units of (seconds)^/^; aj is dimensionless for a normalized function / which has 
dimension of sec^/^. 

Viewing a single-mode of a path as a transmission line, we need an operator corresponding 
to voltage (analogous to the electric-field operator of quantum electrodynamics). 

3. Quantum mechanics stripped of space and time 

Often it is convenient to make a preliminary analysis that skips all the integrals over fre- 
quency by treating quantum states in a toy Hilbert space of finite dimension, which means that 
space and time are collapsed to zero dimensions. (That still leaves polarization, for example.) 
This procedure is equivalent to an analysis allowing for frequency for certain calculations, 
namely when the frequency functions involved are all mutually orthogonal. An example is Eq. 

APPENDIX B: OPERATOR LEMMAS 

For any two operators A and B let [A, B] =^ AB — BA. 
Lemma: For any operators A, B,C, 

[A,BC] = [A,B]C + B[A,Cl 

[AB,C] = [A,C]B + A[B,C]. (Bl) 



60 



Lemma: If [B,C]= 0, then 



[[A,B],C] = [[A,C],B], 

[C,[B,A]] = [B,[C,A]]. (B2) 

Lemma: For any four operators A, B,C, D, 

[AB, CD] = A[B, C]D + [A, C]BD + C[A, D]B + CA[B, D]. (B3) 

Lemma: For n = 1,2, 

n 

[A, 5"] = J2 B^'^[A, B]B''-\ (B4) 

k=l 

Lemma: If A\0) = B\0) = 0, then 

{0\ABA^B^O) = {0\A[B,A^]B^\0) + {0\[A, A^][B , B^]\0) . (B5) 
Lemma: For any two operators A and B\ if [[[A, B''], B''],B'^] = 0, then 

[A, B^''] = n I !lZ:l5t(n-2) ^^^^ ^t] ^ ^t] + QKn-l) [^^ ^t] | ^ (gg) 

Proof: [A, B^'"] = Yl=i B^^'"^^ [A, 5t]5t("~^) and 

[A, = + [[A5t],fit(— fc)] 

whence the lemma follows. □ 

Lemma: If [[[A, A^,A^,A^ = and A|0) = 0, then 

(B7) 

Proo/: Notice that (0|A"At"|0) = (0|A"'i[A, At"]|0) and use Lemma (lB6ll. □ 
Lemma (IB7I) shows how repeated commutators work their way into (0|y4"y4'''"|0). 
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Lemma: Given A\0) = B\0) = and [[A, B% B^ = [A, [A, B^] = 0, then 

(0|A"5^"iO) = n(0|v4"-^fi^("-i)[A,5^]|0), (B8) 
from which it follows that: 

Lemma: Given A\0) = B\0) = and [[A, B% B^ = [A, [A, B^] = 0, then 

(0|[A",S^"]|0) =n!(0|[A,S^]"|0). (B9) 

Note: if A = ttf, then n\~^^'^A^^\0) is an n-photon state; thus two n-photon states of this type 
have as their inner product the n-th power of the inner product of the corresponding 1 -photon 
states. It follows that a unitary transform can convert an n-photon state into a tensor product of 
an (n — l)-photon state and a 1-photon state. 

Lemma: For any two operators A and B such that [B, [A, B]] = and = 0, 1, 2, . . . , we 
have 

[A, B"] = n[A, B]B"-\ (BIO) 

(Proof follows by induction, using Lemma (IBll) .') 
Lemma: For any two operators A and B such that [B, [A, B]] = 0, we have 

[A,exp{B)] = [A,B]exp{B). (Bll) 

(Proof by expansion of exponential, using Lemma (IB3I) .) 
Lemma: Given any operators A, Bi, . . . , Bn, 



A,l[B, 

. i=i 



n /e-1 \ / n 

Y.[IIbA[A,B,]1 Ub,], (B12) 

e=i \j=i J \j=e+i 



with the convention that for any m > n and any Xj 

n 

l[X, = l (B13) 
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1. Implications of commutation rules 

Assume for the rest of this appendix the commutation rules Eqs. (12.61) and (I2.7I) . Then we 
have 

Lemma: For a set of frequencies to, cui, . . . , tUn, with n > 1, 



a{Lj),Y[a\uj) 
i=i 



^5(cu-a;,) Ha^K). (B14) 



j=l k=l 
k^3 



Lemma: Let Sn be the permutation group on 1, . . . , n; then 

m n n 

j = l k = l TV&Sn j = l 

From this follows 
Lemma: 

(0|(fif^:a"')(/i„:a'f")|0) = ^^^n! J duji ■ ■ ■ dUn gl,{uJi, ■ ■ .,uJn)S{uJi, . . .,ujn)K{uJi, • • • ,u;n), 

(B16) 

where S is defined in ( I2.20I ). From this follows another useful fact of norms: 
Lemma: If hn{uji, ■ ■ ■ , tOn) is symmetric under all permutations of its arguments, then 

(0|(/i::a")(/i„:a^")|0) =n! ^ d^i ■ ■ ■ . . . , ^^T- (B17) 

We also have the following relation that allows the calculation of some detection probabilities: 

Lemma: If [b{uj), a^{uj')] = and hn{uji, ■ ■ ■ , tOn) is symmetric under all permutations of its 
arguments, then 

{0\{h::a'b--'){hn:a^'b^^"-'^)\0) = ^^^^^^ (0|(/i::a")(/i„:at")|0). (B18) 

Proof: 

(0|(/i::a'=6"-'=)(/i„:a^'=6^("-^))|0) 



dui - ■■ dUn duj[-- ■ dJ^ hl{uJi, . . . , u;„)/i„(cj^, . . .,uj'J 
{0\a{uji) ■ ■ ■ a{ujk)b{ujk+i) ■ ■ ■ b{ujn)a\uj[) ■ ■ ■ a\ujip{ui^^) ■ ■ ■ 6^^010)- 
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Because a and b work on distinct tensor-product factors, we have 

{0\a{ui) ■ ■ ■ a{uJk)b{uk+,) ■ ■ ■ b{uJn)a\u[) ■ ■ ■ a\uj',)b\u',^,) ■ ■ ■ b\u'J\0) 

= (0|a(a;i) ■ ■ ■ aMa\u[) ■ ■ ■ at(u;;)|0)(0|6(^,+i) ■ • • bMb\u',^,) ■ ■ ■ b\u'J\0). 



The lemma then follows from the symmetry of /i„ together with Lemma (IB15I) . □ 



Concerning broad-band coherent states, from Lemma (IB 141) follows: 
Lemma: For n > 1, 

aH(4)"iO) = nf{u){a\r-'\0). (B19) 
Lemma: For the coherent state defined by Eq. (12.291 ). 



ag\a,af) = J duj g*{iu)a{uj)\a, aj) = a y^J dcu g*{uj)f{iu) j \a,af). (B20) 
Lemma: ForP„(ai,a2 ) defined in Eq. (IT8b . 

oo „ 2 

^n_P„(ai, 02) = I duo^^a}-{ijj)aj{(jj). (B21) 

n=0 j=l 

Proof: By the definition of Eq. ( 13.81) . we have 

00 00 n. 

^nP„(ai,a2) = ^ ^ Pfc(ai)^'n-fc(a2) 

n=0 n=0 k=0 

00/00 00 \ 

= ^ kPk{a,) J2 Pn-k{a2) + Pfe(ai) ^^(n - A;)P„_fc(a2) J . (B22) 

fc=0 \ n=k n=k / 

The lemma then follows from Eqs. (13.41) and (13.111 ). 

APPENDIX C: ALGEBRA OF FREQUENCY-ENTANGLED OPERATORS 

We want to evaluate expressions of the form (0|Pol^Pol|0), where Pol is a polynomial in 
annihilation operators. The general method of evaluation is to use commutation relations to 
rearrange the operators so that, in the end, nothing is left but a number. The commutation 
relations amount to an algebra, which we now construct for the simplest quantum models 
that show how polarization entanglement combines with frequency entanglement. The models 
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cover bi-SU(2) invariant states built up from polynomials in operators of the form g'.Ojbl. 
Consider some number Na of modes derived from a{u;) and of modes derived from b{u;): 

aj(u;),bk(u;), for j = 1, . . . , A^„, and /c = 1, . . . , A^b, (CI) 

with 

a,H|0) = bk{u)\0) = {0\a]{u) = {0\bl{u) = 0. (C2) 
The commutation relations are 

[aj{u;),al{u;')] = 5jk5{u; - u;'), 
[bj{uj),bl{u')] = 5jkS{(jJ-cv'), 

[aj{u), bliuj')] = K(cu), akiuj')] = [bj{uj), bk{uj')] = 0. (C3) 



1. Arrow notation for frequency dependence 

Let p and q denote any of these improper creation or annihilation operators. We consider 
operators of the form 

p.^.q^ J j du duj p{uj)g{uj,uj)q{uj). (C4) 

Let and 6^ denote any of the improper creation operators and fix a square-integral function 

51(0;, a; ); we are interested in the commutator algebra generated by 

a) ■■^■h^ (C5) 

and its adjoint, which is 

(at . ^ . 6t)t ^i,.C -a^a-i -b. (C6) 

We now develop this arrow notation. For any function of two variables g{xi, X2) that en- 
ters as a factor in an integrand, we write 5^ as if the second variable is identified with a 

variable in a following factor or if the first variable is identified with a variable in a preced- 

g 

ing factor; we write <— if the first variable is identified with a variable in a following factor 
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or if the second variable is identified with a variable in a preceding factor. This makes it 
easy to express compound convolutions that will occur in commutators, such as X2) — 

j Jdx dx' f{xi,x)g{x', x)h{x', X2) and a;2) = j J dxdx' f{xi, x)g{x, x')h{x', X2); be- 

cause the order is different in the middle factor, Gi and Gu are distinct. Writing a "•" for 
integration, we then diagram 

^ / N f 9 h 

Gi{xi,X2) as xi ^ ■ ■ ^ X2, 

Gu{xi,X2) as xi X. • • A (CI) 

Altogether there are eight such functions, corresponding to the eight ways to orient a sequence 
of three arrows. 

A compound function such as Gi can itself enter another compound. If G is defined by 

G 

an arrow diagram with the above procedure, then is obtained immediately as that diagram 

G 

while <— is obtained by left-right reflection of the diagram, as in 



We abbreviate repetitive patterns by an exponent; for example 

9 h g h g fg h\'^ g 



(C8) 



(^•A)-^. (C9) 



We define an operation 'Loop' that produces a number from a string of arrows by joining the 
two terminal points; for example 

Loop^-^ • = (^*,/i)= J J dcu dw g{u!,cu)h{uj,cv). (CIO) 

Lemma: 

A.(^.A)" ^ = (cii) 

Lemma: With these arrow rules, if [p, q] — 0, then 

p-^-q^q-^-p. (C12) 
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2. Algebra rules 



Equations (IC3I) imply a commutator algebra generated by 

^■bk and al-^-bl (C13) 



The commutator of these, along with commutators of commutators, etc., generate new fre- 
quency functions. Regardless of the frequency function g, the operator aj ■ ■ b^ acting to the 
right on the vacuum state annihilates the state. We will speak of such an operator, regardless 
of its frequency function, as being of annihilation type. In addition to the creation and annihi- 
lation types that we start with, the commutation relations engender two more types, one type 
has the form aj ■ — ■ or fej ■ ^ ■ bk, the other type is just a number. The point is to evaluate 
expressions of the form (0|Pol^Pol|0) by using the commutator algebra to transform this to a 
form (0|x|0), where x is just a number, extracted in the last step from the normalization relation 
(0|0) = 1. 

The commutation relations among all these types are defined by the following and their 
hermitian conjugates: 



^ ■ Ofc, a; ■ ^ ■ O^J = Ofem aj, ■ ^ ■ ^ ■ aj + djib^^ ■ *~ ■ ^ ■ 6^ 

+ Mfc„Loop(A-^), (C14) 

: 5,,a]-^-^-bl, (C15) 

-■ Okiay ■ ^ ■ am - Oj.m a'l, ■ ^ ■ ^ ■ ak, (C16) 

^ 5kmal-^-^-b], (C17) 

'■ Skiby ^ ■ ^ ■ bm - Sjmbl ■ ^ ■ ^ ■ bk. (C18) 

These relations hold for any integrable functions g{uj, lu') and h{uj, uj'). 





9^ 


■ ttk, a\ ■ 


A 


■bl 


a] ■ 


9^ 


■ flfe, o.£ • 


A 






9^ 


■ bk, a\ ■ 


A 


■ 6t 

m 




_ 9^ 


■ bk, b\ ■ 


A 


■b^ 

m 
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3. Subalgebra for fixed g 



For constructing examples, we will use a subalgebra in which a chosen g(uj,uj') plays a 
distinguished role. Given any such g, define for n = 1, 2, . . . , 



9 9*V^ 9 



9_,(^_9_ 

^ . Cy 

9 



n-l 



9n^ 

gn 

By writing out a few arrow expressions, one shows for positive integers r and s 



9s ^ . ^ 
,9s _ ^ 

9s ^ 

9s ^ . £l 

Gn 



9^ . £y+*-l ^^ r+s- l 



9r+s 
9r+s 
Gr+s 

n 



With these and the definition in Eq. (IC19I) . we find 



9; 



Gr 



G'f 



■ bk, ■ 



t 9s 



t 9s 



G. 



■hi 



■hi 



a„ 



c- t Gr+s-l e- ,+ 

0km O-i ■ * ■ O-j + Oji 0^ ■ 

+ 5je6km^OOp{Gr+s-l), 



Gr+s-1 



■ hk 



Ski a'- 



t 9r+s 



hi 



Gr 



Gr 



Oke (^i ■ ■ O-m ~ Ojm Q-e ' ^ 



0^ ■ > ■ Ok, a. ■ —>■ ■ 0' — Okm 0,1, ■ — > ■ 0. 



(C19) 



(C20) 
(C21) 

(C22) 
(C23) 
(C24) 
(C25) 
(C26) 
(C27) 



'km Of 



(C28) 
(C29) 

(C30) 

(C31) 
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bk, b\ ■ ■ b^, = 5keb\ ■ 



Gr+S T _ r , t 



■ bk 



(C32) 



4. Commutators of powers of operators 



To evaluate powers of operators, we shall need some repeated commutators. 
Lemma: For any positive integers r,s,q. 



■ bk, a\ ■ ■ bl] , al ■ ■ bl] = 26ip6km a\ ■ ■^'^t!±| ^ -bL- 



(C33) 



It is worth noticing that this double commutator is of creation type, and thus commutes with 
all other operators of that type. The commutation relations also imply 

Lemma: Regardless of what functions decorate the arrows, 

[[a • ^ ■ bj, ■ ■ bll ■^■bl]= 0, if j ^ k. (C34) 

In particular the operators ■ ^ ■ bk and a] • — ^ ■ above fulfill the conditions for this 
lemma. 

Proposition: Write {abj) for a ■ ■ bj and for non-negative integers ni,n2, . . . ,nK, let 
n = Y^j^iUj-, then 

n(«^^)"^) (n^""^^-)^"^) = {Q\{ab,r{ab,r\Q). (C35) 

Proof: The proof employs the symmetry operator of Eq. ( 12.201) . The left side of the equation 
is the inner product of a vector |0) with itself, where 

10) = /i2n:Pol^|0) 

)Pol^|0), 



with 



and 



ni 



poit= (n«^(^.w( 



(C36) 



(C37) 
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Because of the commutativity relations, the integral is unchanged under certain permutations 
of the arguments of /i2n- Interchanging the operators a^{ujj) among themselves yields 

10) = {S{uJi,...,UJn)hn)M^O). (C38) 

In addition, letting Cuj = {luj, luj), we have 

\<P) = {S{u^,...,uJ^)h2n)■M^0). (C39) 

These two symmetries for and its product form defined in Eq. (IC36I) imply a third symme- 
try: 

|0) = (5(cDi,...,cD„)/i2n):Polt|O). (C40) 

From this symmetry and Lemma ( IB15I ) the proposition follows. □ 
This generalizes Lemma (IB18I ). 

We shall need to refer to a function 



E^in)"^ - {0\{abnaby-\0). (C41) 

This function is independent of the choice of operators a{uj) and b{uj) as long as each oper- 
ator satisfies the commutation relations Eqs. ( I2.6I ). ( I2.7I) and the two operators are mutually 
orthogonal. 

q* 

Proposition: Write (ah) for a ■ — > ■ 6. Then 

= ^(0|(a6)"(a6)t"|0) 
nl 

= V ■ -Y\\Loop(G.j)Y\ (C42) 

ui,...,u„=0 j=l 
Y, ki/k=n 

where the Uk are restricted as indicated. 



1 



Proof: From symmetry considerations and Lemma (IB 151) we get 

1 1 /" _ 

{Q\{ab)"'{aby"'\0) = — I duidCbi - ■ ■ dujndunduj[du}[ - ■ ■ duj'^du'^ 



[g{uJi, UJi)--- g{uJn, UJn)] ^[) ■ ■ ■ 9*{^'n^ ^'n)] 



X 



VTreSn j = l / \lT&Sn j = l 
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By virtue of Eq. ( IC39I ). any permutation of the variables uj^^^ is compensated out by a corre- 
sponding permutation of the variables lo'^^, which implies 
1 



= n! y duidi^i - ■ ■ dunduJn\S{uJi, . . . ,UJn)g{uJi,UJi) ■ ■ ■ g{uJn,UJn)\ ■ 

(C44) 

The effect of each of these remaining permutations is to generate a product of integrals, each of 
convolutions of (jf's with matching convolutions of g*'%, according to the cycle structure of the 
permutation group Sn- The convolutions generated by a permutation are those that correspond 
to the cycles of its conjugacy class. Each conjugacy class is characterized by some [z/i, . . . , z/„] 
where vi is the number of one-cycles, U2 is the number of two-cycles, etc. (with X]fe=i = n) 
The number of permutations in a conjugacy class is just that stated in the Proposition 

□ 

Examples are 

2,(0) = 1, 

2,(1) = Loop(Gi), 

S,(2) = [Loop(Gi)]2 + Loop(G2), 

S,(3) = [Loop(Gi)]=^ + 3[Loop(Gi)][Loop(G2)] + 2Loop(G3). (C45) 

As a check on the Proposition, these examples can also be demonstrated by repeated use of Eq. 
(IB7I) along with the commutation relations Eqs. (IC28I) - (IC32I) . 

We will need to deal with partial traces, as defined in Sec. El in particular we need 



Tr„[(a6)^'^|0)(0|(a6)"] = dCjdCj' 



dojdoj'giu:, Cj)g*{oj', lj') (0„|a"(a;')at"(a;) |0„) 



n\ / dujdu)' 



j du; g{u;, u;)S{u:)g*{u;, 



6t"(cD)|05)(05|6"(Lj') 



6t"(c:;)|0b)(0,|6"(cD'), (C46) 
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where we have used the compact notation defined in Sec. El with g{(jJ, uj) =^ g{oJi, ui) ■ ■ ■ 
g(uJn,uJn), and the second equality invokes Lemma (B15). Because of symmetry of the b 
operators under permutations of their arguments, this simplifies to 



TrJ(a6)1-"|0)(0|(a6)"] I dC^dC^ 



dujg{uj,uj)g*{uj,uj') 



6t"(i,)|0,)(0,|6"(cD'). 



(C47) 



Writing out the compact notion, we note that the inner integral above is 

jdu}g{uj,uj)g*{uj,uj') = Yl(^jduJ g{uj,ujj)g*{u,Uj] 



n 



g g ~, 

< >LJ^ 



(C48) 



As a check, we note that the trace over the 6-mode of Eq. (IC47I) results in Proposition (C42), 
as it should. 



5. Examples of frequency functions 

We consider a family of functions g(;{uj, u) and show two limiting cases. For any real-valued 
functions (p{uj) and (p{uj) and positive real parameters a and a, let 



'(TO" 



UJ — UJq LU — UJq 



a 



a 



(C49) 



(C50) 



where we define 

^ ' '^/e + ~^ + c){x+yf + (ye + i + c)~\x-y) 

= y|exp|- + 1 (x' + y^) + 2Cx2/] }; 
regardless of the value of C, 

/oo 
dxdy\F{C;x,y)\^ = l. (C51) 
'OO 

Thus for any choice of center frequencies c^o and ujq, bandwidth parameters a and a, and phase 
functions 0(u;) and 4>{Cj), we get a family of (7^'s. 
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By changing integration variables, one shows 



Loop 



9<, 9l 



Loop 



F F' 



(C52) 



(C53) 



and similarly for the convolution of (7^'s 

It remains to compute the convolution integrals for F. Although gi^(uj,uj) need be neither 
real-valued nor symmetric in to and cu, the function F{(; x, y) is real and symmetric, implying 

= (C54) 



For this reason arrow expressions built up from convolutions of factors of F are invariant under 
any number of reverses of arrow directions; to emphasize this indifference to arrow direction 
in F (but not (7^), we write 

. (C55) 



To compute the convolution integrals for F, consider a sequence of Q. Abbreviating F{C 



by Fj and y + 1 by cj, we have 

X < — > ■ < — > y 

def 



- [ dz exp{-[c,(a;^ + z'^) + Ck{z'^ + y'^) + 2CjXz + 2C,kzy]} 

TT J 



Cj + Cfc V VT 



exp < - 



Cj + Cfc 



2 -1 1/2 



x' + y') - 



Cj + Cfc 



xy 



2 / CiCk 

F\ \x,y 



Cj -\- Ck \ Cj -\- Ck 

Setting y = X and integrating yield 



(C56) 



F- F- 
(Vj) Loopf^-^ 



1. 



(CSV) 



From this and Eqs. (IC50I) and (IC56I) follows the 
Lemma: 



a^l'SSj y = ^yK{C,n) F{Cn;x,y), 



(C58) 
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where we have written F(^) as an abbreviation for F{(;-,-) and 



Ci = c, 

Cn+1 = — 



ClCn 



2 



From the lemma and these equations we find 



Loop ^— ^ 



Loop 
Loop 
Loop 



(m no 

/F(c)y (F{o 



with K{C,n) as defined in (IC62I ). Similarly this and Eq. (IC19I) yield 



Loop 



Note that Eq. (IC62I) implies 



(Vn>2) ft:(C,n)<— . 



(C59) 
(C60) 
(C61) 

(C62) 



(C63) 



(C64) 



(C65) 



For more efficient calculation, define n) for n > 2 by 

K(C,n) = 4^£ii. (C66) 

From Eqs. (IC59I) and (IC62I) . it follows that i?(C, n) is a rational function of C satisfying the 
following recursion relation 



Lemma: 



i?(C,2) = L 
(Vn>2) i?(C,n+l) = ( 1- 



4(e + 1; 



(C67) 
(C68) 
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From Lemma (IC58I) and Eq. (IC53I) one also finds 

— ^ = k2n+i —1=^ e '^^ ^ F ( C2n+i; 
veer 



a;i— = cr ^A;2„ exp{i[0(u;i) - ^(cua)]} F ( C; 



0" cr 



(7 



a 



(C69) 
(C70) 
(C71) 



6. Limiting cases 



Case I. No frequency entanglement: giiuj.u)) = f(uj)h(Co), normalized so that 
J J dujdCu \g{uj,Cj)\'^ = J duj 1/(^0") ^ = J duj \h{u)\'^ = 1. In this case, one skips the fancy 
commutation relations because the operators all factor; we find 

Lemma: 

5i(n) 1 (0|(a ■ ^ ■ 6)"(at ■ ^ ■ 6t)"|o) = n\, (C72) 

which follows from the discussion of broad-band coherent states in Sec. 12 CI 
Case II. Limit as ^ — > ±00. It makes no sense to ask for the limit of (?^; however, we explore 
large values of \(^\ by looking at the limit of the commutation relations. From Eqs. (IC62I) and 
(IC64I) we see for this limit 

(forn > 2) lim fc„ = lim Loopf— = lim Loopf^^ • 

resulting in specializing the commutation Eqs. (IC28I) - (IC32I) . for sufficiently large \(^\, to 

3* h t ^ At 



(C73) 



lim 



t 9 



bl 



lim 

C^itoo 



lim 
lim 



Gi f Gi 
b] ■ ^ ■bk,al - - bl 



Gi 



bk, bl 



Gi 



0, 
0, 
0, 
0. 



(C74) 
(C75) 
(C76) 
(C77) 
(C78) 
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is effectively 



In this case, the double commutator ■ ^ ■ 6^, ■ ^ ■ &m 5 '^1 ' ^ ' 
zero, so that Lemma (IB6I ) applies. 

From Lemma ( IB9I) and Eqs. ( IC74I )-( IC78I ) follows the corresponding rule for evaluating 
Case-II operator products: 

Lemma: 

Eu(n)= lim ^(0\(a-C-by(a^ ■-^■bA"\0) = l. (C79) 
C^itoo nl \ /V / 

APPENDIX D: FOURIER TRANSFORMS IN SPACE AND TIME 

Let f{x, t) be any operator- valued function for which Fourier transforms make sense, and 
define the Fourier transform pair: 

/oo 1*00 
dt / c/x/(x,t)e*("*+'=^\ (Dl) 
CXD J —00 

/OO poo 
duj / rfA;7(cu,A;)e-'("*+'="). (D2) 
-00 J —00 

On taking the adjoint of these equations, one sees that: 

Lemma: The Fourier transform of P{x, t) is related to the adjoint of the transform of /(x, t) 
by 

f{-uj,-k) = U{uj,k))\ (D3) 

Lemma: If /(x, t) = p{x, t), then 

J{-u;,-k) = (J{u;,k))\ (D4) 

so that the Fourier transform of a hermitian operator function is specified for the whole (to, k)- 
plane once it is specified for any half -plane touching the origin. In particular Eq. (ID2I) can be 
replaced by 

/oo 1*00 ^ 

duj dA; (^7(u;,fc)e-'("*+'=")+7^(cj,A;)e*("*+^")j , (D5) 

where the integrand is integrated over the half-plane A; > 0; alternatively f{x,t) can be ex- 
pressed by the same integrand integrated over the half -plane defined by u > 0. More generally, 
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the region of integration needs to be obtained from the region excluded by reflection through 
the origin; thus the region of integration need not be aligned with the axes and need not have a 
straight boundary. 

Lemma: If f{iu, k) = f\uj, k), then 

J{-u,-k) = (j{u:,k))\ (D6) 
Now consider the form of a Fourier transform of any solution to the wave equation 

{dl-dl)g{x,t) = Q- (D7) 
this equation has as its general solution 

g{x,t) = g+{t-x)+g+{t + x), (D8) 
so that the transform of g has the form 

g{uj, k) = [5{k + uj)g^{u) + 5{k - uj)g_{uo)l (D9) 

where 

g^[uj)^^-= due'^-g^iu). (DIO) 



oo 



/27r 

Taking the inverse transform of Eq. (ID9I) . one obtains for the general solution to Eq. (ID7I) 

g{x,t) = / du;[g4cu)e-'-('-^^ +g_iu;)e-'^^'^% (Dll) 

V 2lT J -oo 



APPENDIX E: EXPANSION OF LIGHT STATES IN TENSOR PRODUCTS OF BROAD- 
BAND COHERENT STATES 

Consider the subspace of broad-band coherent states spanned by aj"|0), n = 0, 1, 2, It 

follows from Louisell O, P- 106] that the unit operator on this subspace is 

\a,af){a,af\ . (El) 

TT 

Next, let F be any set of orthonormal functions fj{uj), j = 1,2, This implies that the set 

of operators \aj, af.){aj, af .\, j = 1,2, ... , are mutually orthogonal projections. We say a set 

77 



of light states is 'coherently expressible with respect to F' if for some set F of orthonormal 
functions fj{uj),j = 1,2, . . . , every state of the set is some weighted sum (or integral) over j 



and aj of states of the form 



ni«„4,), (E2) 



where the product is a tensor product. The unit operator for the vector space of such states 
coherently expressible with respect to F is then 

J \aj,af^){a,af^\——. (E3) 



Subtleties of coherent states in infinite dimensional vector spaces are touched on in 1118 
503-512]. 



pp. 



APPENDIX F: MATLAB PROGRAMS FOR SECTION 10 

Here we record the MATLAB scripts and functions used to generate Fig.|71 

(1) Partfn.m starts the calculation by preparing a list yList of the partitions of integers 
needed for Eq. (IC42I) . It stores yList inafile Part .mat . It needs to be run only once, with 
a value of nmax = floor (3 * mujnax + 16) , where mu_max is the highest value 
of jj, covered. (For nmax = 32, Partfn.m takes 20 minutes on a Pentium-4 desktop 
computer.) 

(2) muRun . m generates data in the cell variable muResult for later plotting. It has a section 
that is easily modified in order to zoom in on one or another parameter region. This part 
takes values for pdark, Vda and r/trans for each detector. It takes the parameter vsq = |f p 
that is the fraction of energy tapped by the eavesdropper. Finally it takes a parameter for 
the order of Renyi entropy considered in calculating the eavesdropper's entropy. (Caution: 
stronger eavesdropping attacks are expected to be implemented in the future.) In order to 
speed calculations, it calls Xif n . m to pre-compute values ofE((^,n), and it also computes a 
list of all n\ for n = 0, . . . , nmax . 

(3) EntPlts.m plots families of curves, such as that shown in Fig. |71 working from 
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muResult , obtained either from immediately prior running of muRun . m or from loading a 
previously saved muResult . 

The rest of the scripts and functions are called directly or indirectly by muRun . m : 

(4) muScript . m is a macro called three times by muRun . m . 

(5) Xif n .m pre-computes a list of values of n) and stores them in XiAr , in order to 
speed the computation of probabilities. 

(6) AvEntf n .m computes Eq. (IIO.IOI) . 

(7) Tf n .m supports either (110.591) or (110.601) by computing (110.391) or (110.551) . respectively, 
as specified by CaseCode , using Ff n . m . (So far, the only energy distribution implemented 
is the Poisson energy distribution, for which C„ = e~^n'^/n\.) 

(8) Ff n . m computes either (110.491 ) or (110. 551 ). as specified by CaseCode . 

(9) Gkmfn.m is called by Ffn.m to compute Eq. (110.491 ): it uses (global) XiAr, 
CaseSetUp , and gamTab set up by muRun . m . It calls sumfn.m. 

(10) sumf n .m is a summing routine for efficient calculation of sums in which the ratio of 
successive terms is pre-computed. 

(11) Gnf n . m is called by Ff n . m to compute Eq. (110.581) . 

(12) Gmu_kmf n . m computes Eq. (110.611) for the special case of C„ = e~^/i"/n!. 

(13) Gmuf n . m computes Eq. (110.631) for the special case of C„ = e~^/i"/n!. 

(14) GCf n . m. Not yet implemented. 

Here is the code for all these except GCf n . m, not yet implemented. 
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MATLAB PROGRAMS 

(1) Partfn.m 

function [ ] = Partfn(nmax) 
% 26 OCT 04 

% sets up cell for partitions from n—lton — nmax 
% assumes nmax > 2 

% test and run for nmax = 3 2 on 26 Oct 04 (took 20 minutes) 

yList = cell ( 1 , nmax) ; 
yList{l} = [1] ; 
yList{2} = [0 1; 2 0] ; 
for jcell = 3:nmax 

yList{jcell} = Upf n (yList{ jcell-l}) ; 
end 

save Part. mat yList % ESSENTIAL RESOURCE! 

function [y] = Upfn (Ar) 

% Used in generating partitions of n + 1 from partitions of n; 
% Ar has a row for each u vector of partitions of n. 

temp = size (Ar) ; 

oneCol = ones (temp ( 1 ), 1 ) ; 

zeroPad = zeros ( size (Ar) ) ; 

oneColPadded = [oneCol zeroPad] ; % row width is n + 1 

zeroCol = zeros (temp ( 1 ), 1 ) ; 

ArPadded = [Ar zeroCol]; % add a column to get n + 1 columns 

ArPl = ArPadded + oneColPadded; 
% Up proper 

for jRow = l:temp(l) 

for j = 1 : temp (2 ) 
if Ar ( jRow, j ) > 

Vec = ArPadded ( jRow, :) ; 
Vec ( j) = Vec ( j) -1; 
Vec(j+1) = Vec(j+1)+1; 
ArPl = [Vec;ArPl] ; 
end % of If 
end % of For j 
end % of For jRow 
ArPl = sort rows (ArP 1 ) ; 
tempz = size (ArPl); 
y = ArPl (1, : ) ; 
for jRow = 2 : tempz (1) 
if ArPl ( jRow-1, : ) == ArPl (jRow, :) 
else 

y = [y;ArPl (jRow, : ) ] ; 

end % of If 

end % of For jRow 
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(2) muRun . m 

% Compute and store and plot /^-dependence of figure of merit 

% and pGood, pSiftErr, and AvEnt 

% Define figure of merit = pGood . * AvEnt 

% 1 NOV 04 

global alphaO dark vsq CaseSetUp zetaVec XiAr NsqFac yList 
^ ********************** VARY TO DEFINE CASE ********************** 

etaDet = [0.1 0.1 0.1 0.1] ; 

% etaDet =[1111] ; 

etaTrans = [1 1 0.1 0.1]; 

% etaTrans = [1 1 1 1 ] ; 

dark = 5.*10"(-5) .* [1 1 1 1]; 

% dark = [0 ] ; 

alphaO = 1-etaDet .* etaTrans; 

vsq = .25; % vsq is fraction of Bob's energy tapped by Evang. 

% vsq = ; 

zetaVec = [1 10 100 1000 ];% Row vector ( 1 , length ( zetaVec) ) ForCase 

Ry = 1.1;% Order of Renyi entropy 

% can put in knob variables later 

% need to make global array of |C„p if cases 5, 6 used. 

% parameters for range and fineness of 

fine_incr = .0001; 

mu_begin = ; % may add a little to avoid problem when Pdark = 0; 
mu_f ine_max = .007; 

n_fine_max = floor ( (mu_fine_max - mu_begin) /f ine_incr ) ; 
mu_max = 0.04; 

incr = . 02; % increment /i after first steps of incr/n_f ine_max 
n_incr_max = 1 + floor ( (mu_max - mu_begin - 
n_f ine_max*f ine_incr ) /incr) ; 
% mu_max = (n_incr_max+l ) *incr 

muVec = zeros (1, n_incr_max+n_f ine_max+l) ; % will be loaded with values 
^ ************************ gj>^X) OF CASE DEF ************************ 

CaseSetUp = cell (1,6); 
CaseSetUp{l} = etaDet; 
CaseSetUp{2} = etaTrans; 
CaseSetUp{3} = dark; 
CaseSetUp{4} = vsq; 
CaseSetUp{5} = zetaVec; % row vector 

CaseSetUp{6} = Ry; % Renyi entropy of order Ry . load Part .mat 
nmax = 32; 

if ( (zetaVec == zetaVecOld) & (nmax == nmaxOld) ) 

else 

XiAr = Xifn ( zetaVec, nmax, yList ) ; % assumes yList onhand 
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% XiAr = cell (1, nmax+1) ; 

% cell{n+l} is Xi (zeta, n) , a column vector (1, length (zetaVec) ) 
zetaVecOld = zetaVec; 
nmaxOld = nmax; 

end 

% set up table of factorials to speed calculation 

global gamTab 

gamTab = cell ( 1 , nmax+1 ) ; 

for kk = 1:33 

gamTab{kk} = factorial (kk-1 ) ; 
end 

% END of setup. m 

% GET LIMITS for muMeritLims 
CaseRestore = CaseSetUp{5}; 
CaseSetUp{5} = 0; 
muScript 

pGoodZO = pGood; 
AvEntZO = AvEnt; 
FigMerZO = FigMer; 
pSiftErrZO = pSiftErr; 
CaseSetUp{5} = 10'^50; 
muScript 

pGoodZInf = pGood; 
AvEnt ZInf = AvEnt; 
FigMerZInf = FigMer; 
pSiftErrZInf = pSiftErr; 

% Get in-between values of zeta 
CaseSetUp{5} = CaseRestore; 
muScript 

muResult = cell (1,5); 
muResultjl} = CaseSetUp; 

muResult{2} = [pGoodZO;pGood;pGoodZInf ;muVec] ; 
muResult{3} = [AvEntZO; AvEnt ; AvEnt ZInf ;muVec] ; 
muResult{4} = [FigMerZO ; FigMer ; FigMerZInf ; muVec] ; 
muResult{5} = [pSiftErrZO;pSiftErr;pSiftErrZInf ;muVec] ; 

% next can be run separately as Ent P 1 1 s 

FigMerPlt = muResult{4}; 
plot (muVec, FigMerPlt (2, : ) ) 

xlabel ( ' mu ' ) 

ylabel ( 'FigMerit ' ) 

t it le ( ' FigMerit vs. mu; zeta 

text (.5,. 05, '\itp_{\rm dark} 
det} = . 1, \eta_{\rm Trans} =1.') 



= ; zeta = inf -- ' ) 
= 0, \eta_{\rm 
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(3) EntPlts .m 

% Program For plotting data In mu Re suit gotten from running muRun or 

% loading muResult 

% Plots choice of "kind". 

% 6 NOVEMBER 2004 

kind = input (' 1 for pGood, 2 for AvEnt, 3 for FigMer, 4 f 
pSiftErr ' ) 

% takes iJ, values as vector muX, an edited muVec taken from edited muResult 
pltVec = muResult{kind+l}; 

% ******************* VARY UNTIL next asterisks TO EDIT 

startDim = size (pItVec) ; 

% pItVec ( : , 4 : startDim ( 2 ) ) = [];%** TEMPORARY; comment OUT 

pItVec ( : , 1 :20) = [ ] ; 

^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

tempDim = size (pItVec) ; 

mux = pItVec (tempDim(l) , :) ; % retrieves (edited) muVec 
% kind = 1 gets pltVec = pGood 
% kind = 2 gets pltVec = AvEnt 
% kind = 3 gets pltVec = FigMer 
% kind = 4 gets pltVec = pSiftErr 

CasePlt = muResult{l}; % CaseSetUp For data In muResult 

etaDetPlt = CasePlt{l}; 
etaTransPlt = CasePlt{2}; 
darkPlt = CasePlt{3}; 
vsqPlt = CasePlt{4}; 
RyPlt = CasePlt{6}; 

pdarkStr = 

[ ' [ ' , num2str (darkPlt (1) ) , ' , ' , num2str (darkPlt (2) ),',',.. . 

num2str (darkPlt (3) ) , ' , ' , num2str (darkPlt (4) ) , ' ] ' 
etaDetStr = 

[ ' [ ' , num2str (etaDetPlt (1) ) , ' , ' , num2str (etaDetPlt (2) ) , ' , ' , 
num2str (etaDetPlt (3) ) , ' , ' , num2str (etaDetPlt (4) ) 
etaTransStr = 

[ ' [ ' , num2str (etaTransPlt (1) ) , ' , ' , num2str (etaTransPlt (2) ) , 
' , ' , num2str (etaTransPlt (3) ) , ' , ' , num2str (etaTransPlt (4) ) , ' 

Titles = cell (1,4); 

Titlesjl} = 'Prob. of correct, sifted bit vs. mu ' ; 
% factor of 1/2 For prob of correct basis 

Titles{2} = . . . 
[' Evangeline '' s AvEnt for correct, sifted bits vs. mu; R 
' , . . . num2str (RyPlt) ] ; 
Titles{3} = 'FigMerit vs. mu ' ; 

Titles{4} = 'prob. of error in sifted bits vs. mu ' ; 
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Ylab = cell (1,4); 

Ylab{l} = 'pGood'; 

Ylab{2} = 'Evangeline AvEnt ' ; 

Ylab{3} = 'FigMerit'; 

Ylab{4} = 'pSiftErr'; 

% plot (muX,pltVec) 

plot (muX,pltVec (1, : ) ,muX,pltVec (2, : ) , ' — ' ,muX,pltVec (3, : ) , '- 

muX,pltVec (4, : ) , '- ' , muX, pltVec ( 5, : ) , ' - ' , muX, pltVec ( 6, : ) , '- ' ) 
xlabel ( 'mu' ) 
ylabel (Ylab{kind}) 

title (Titles{kind}) 
axis([0 max (muX) -inf inf ] ) 
gtext ( [ ' \itp_{\rm dark} = ' , pdarkStr, . . . 

'\it, \eta_{\rm det} = ',etaDetStr, '\it, \eta_{\rm trans} = 
' . . . , etaTransStr, ' , vsq = ' , num2str (vsqPlt ) ] ) 

% muResult = cell (1,5); 

% muResultjl} = CaseSetUp; 

% CaseSetUp = cell (1,6); 

% CaseSetUpjl} = etaDet; 

% CaseSetUp{2} = etaTrans; 

% CaseSetUp{3} = darkPlt; 

% CaseSetUp{4} = vsq; 

% CaseSetUp{5} = zetaVec; % row vector 

% CaseSetUp{6} = Ry; % Renyi entropy of order Ry 

% muResult{2} = [pGoodZO;pGood;pGoodZInf ] ; 

%muResult{3} = [AvEnt ZO ; AvEnt ; AvEnt ZInf] ; 

%muResult{4} = [FigMerZO ; FigMer; FigMerZInf ] ; 

% muResult{5} = [pSif tErrZO ; pSif tErr ; pSif tErrZInf ] ; 



84 



(4) muScript .m 

% 2 Nov 04 for use by muRun . m 

zVecLen = length (CaseSetUp{5}) ; 
pGood = zeros ( zVecLen, n_incr_max+l+n_f ine_max) ; 
% +1 so ForLoop starts at 1 with /i ^ 

AvEnt = pGood; % same format 
pSiftErr = pGood; % again same format 
% Ry order of Renyi entropy from muRun 
% Do fine steps at beginning 

%mu = .00001;% sloppy fix of problem with /j, = Q when Pdark — 
mu = ; 

for kt = l:n_fine_max 

tOllO = Tfn ( [2 mu] , [0 1 1 0] ) ; 
tlOOl = Tfn([2 mu],[l 1]); 
tlOlO = Tfn([2 mu],[l 1 0]); 
tOlOl = Tfn([2 mu],[0 1 1]); 
pGood(:,kt) = tOllO + tlOOl; 
AvEnt (:,kt) = AvEntfn([2 mu],Ry); 

pSiftErr ( : , kt) = (tlOlO + t 0101 )./ (tOlOl+tOllO+t 1001+t 1010 ) ; 
FigMer = pGood .* AvEnt; 
muVec(kt) = mu; 
mu = mu+incr /n_f ine_max; 
end 

mu = mu- .0001; % don't need fix > 

for kt = n_f ine_max+l : n_incr_max+H-n_f ine_max 
tOllO = Tfn ( [2 mu] , [0 1 1 0] ) ; 
tlOOl = Tfn ( [2 mu] , [1 1] ) ; 
tlOlO = Tfn ([2 mu],[l 1 0]); 
tOlOl = Tfn([2 mu],[0 1 1]); 
pGood(:,kt) = tOllO + tlOOl; 
AvEnt (:,kt) = AvEntfn([2 mu],Ry); 

pSiftErr (:, kt) = (tlOlO + 1 0101 )./ (tOlOl+tOllO+t 1001+t 1010 ) ; 
FigMer = pGood .* AvEnt; 
muVec (kt ) = mu; 
mu = mu+incr; 
end 
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(5) Xifn.m 

function [XiAr] = Xif n ( zeta, nmax, yList ) 

%XiAr(m,n+l) = Xi ( zeta (m) , n) 

% yList in 7matlab/qed/Part.mat 

% (Before running Xif n . load Part .mat) 

% 1 NOV 04 

% 26 October 04 Preliminary test gets OK limits 

% 30 OCT Made XiAr into cell ( 1 , nmax+ 1 ) to get enough dynamic range 

% kappaAr (m, n) = kappa ( zeta (m) , n) 

% Accepts a vector of values of zeta 

% Fails if nmax > nmax used in running partfn.m 

% use in loop for Xi (zeta (m) , 0) ... Xi ( zeta (m) , nmax) by 

%XiAr(m, 1) ... XiAr (m, nmax+1 ) 

zlen = length ( zeta) ; 

% *XiAr = zeros ( zlen, nmax+1 ) ; 

XiAr = ce 1 1 ( nmax+1 ) ; % allows much bigger range of values than array does 

XiAr{l} = ones ( zlen, 1 ) ; 

XiAr{2} = ones ( zlen, 1 ) ; 

kappaAr = kappaf n ( zeta, nmax) ; 

if nmax < 3 

return 

else 

for ncc = 3: nmax+1 
XiAr{ncc} = zeros ( zlen, 1 ) ; 
Lst = yList{ncc-l}; 
dims = size (Lst); 
jmax = dims ( 1 ) ; 
for j = 1 : jmax 

XiAr{ncc} = XiAr{ncc}+termfn (Lst ( j ,:), kappaAr ) ; 

end % For j 
end % For ncc 
end % If nmax 

function [kappaAr] = kappaf n ( zeta, nmax) 
% kappa (zeta (m) , n) 
% 1 NOV 04 

% Accepts a vector of values of zeta 

% kappaAr (m, n) = kappa ( zeta (m) , n) 

% kappa ( zeta, n) = R ( zeta, n) . /sqrt ( zeta''2 + l ) 

X = zeta . " 2 ; 

R = zeros (nmax, length (zeta) ) ; % will be transposed later 
R(2,:) = ones ( 1 , length ( zeta) ) ; 
pvec = X. / (4 . * (x+1) ) ; 

for net = 3 : nmax 

R(nct,:) = 1 . / (1-pvec. *R (nct-1, : ) ) ; 
end 
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kappaAr = R * diag ( 1 . /sqrt (x+1 ) ) ; 
kappaAr(l,:) = ones ( 1 , length ( zeta) ) ; 
kappaAr = kappaAr ' ; % ' 

% checked asymptotic — > 2/ (sqrt (x + 1) +1) as n gets big 

function [yCol ] = termfn (vec, kappaAr ) 

% does Column over zeta (m) For one term where vec is a partition of n 

% For n \le 32 

% using kappaAr (m, n) = kappa ( zeta (m) , n) 

nl = length (vec); 

dims = size (kappaAr) ; 

zlen = dims ( 1 ) ; 

yCol = ones ( zlen, 1 ) ; 

for j = 1 : nl 

yCol = (kappaAr (:, j )./ j ). "vec ( j ) /factorial (vec ( j )). *yCol; 

end 

yCol = yCol . *f actorial (nl ) ; 
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(6) AvEntf n .m 



function [ z_out ] = AvEntf n (CaseCode, R) 
% Assumes CaseCode is [2 mu] or [3 n] 

% 29 October put in "if tlOOl + tOllO > 0" to get rid of 0/0; untested 
% 30 OCT 2004 get zeta from setUp 
global CaseSetUp 
if CaseCode (1) ==2 % BIG BLOCK 
mu = CaseCode (2 ) ; 

kmax = 16+3 . *mu; % from study with poisf n .m 
% numerator = 0; 

zVecLen = length (CaseSetUp{5}) ; 

numerator = zeros ( zVecLen, 1 ) ; 
for k = : kmax 
mmax = kmax - k; 
for m = : mmax 

tlOOl = Tfn ( [5 mu k m] , [1 1] ) ; 
tOllO = Tfn ( [5 mu k m] , [0 1 1 0] ) ; 

if tlOOl+tOllO > 

pEv = tlOOl . / (tlOOl+tOllO) ; 

y = log(pEv."R + ( 1-pEv) . "R) . / ( ( 1-R) *log (2 ) ) ; 
numerator = numerator + (t 1001+t 0110 ) . *y ; 
end%if tlOOl+tOllO . . . 
end % for m 
end % for k 

denom = Tfn ( [2 mu] , [1 1] ) + Tfn ( [2 mu] , [0 1 1 0] ; 
z_out = numerator . /denom; 
else % BIG BLOCK 
% assume [3 n ] 
n = CaseCode (2 ) ; 
numerator = 0; 
for k = : n 
for m = : n-k 

tlOOl = Tfn([6 n k m],[l 1]) 

tOllO = Tfn([6 n k m],[0 1 1 0]) 

if tlOOl + tOllO > 

pEv = tlOOl/ (tlOOl+tOllO) ; 

y = log(pEv^R + ( 1-pEv) ^R) / ( ( 1-R) *log (2 ) ) ; 
numerator = numerator + (tlOOl+tOllO) . *y; 
end%if tlOOl + tOllO... 
end % for m 

end % for k 

denom = Tfn ([3 n],[l 1]) + Tfn ( [3 n] , [0 1 1 0] ) ; 
z_out = numerator . /denom; 
end % BIG BLOCK 
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% function [y_out ] = Tf n (CaseCode, nVec) 
% \mathcal{T} for various cases 

%% nVec is negatively coded bit vector; this restricted implementation fails if nVec has 
%% more than 4 zeros. 

% Casecode can be [1 Cpt] [2 mu] [3 n] [4 Cpt k m] [5 mu k m] 
% [6 n k m] 

% Cpt is a real number or integer that points to an array Cvec of coefficients. 
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(7) Tfn.m 



function [y_out ] = Tf n (CaseCode, nVec) 
% \mathcal{T} for various cases 
% 19 OCT 04 

% Casecode can be [1 Cpt] [2 mu] [3 n] [4 Cpt k m] [5 mu k m] 
% [6 n k m] 

% Cpt is a real number or integer that points to an array Cvec of coefficients. 

% nVec is negatively coded bit vector; this restricted implementation fails if nVec has more 
% than 4 zeros. 

ctab = [0000; 1000; 0100; 1100;... 

0010; 1010; 0110; 1110;... 
0001; 1001; 0101; 1101;... 
0011; 1011; 0111; 1111]; 

% set up Index as vector that shows where the zeros in nVec are located 

% Ct r gets incremented to the number of zeros in nVec. 

temp = size (nVec) ; 

dim_nVec = temp(2); 

dim_Ctr = dim_nVec - sum(nVec); 

j Index = 0; 

if dim_Ctr == % starts big block 

y_out = (-1 ) " sum (nVec) . * Ffn (CaseCode, nVec) ; 
else 

y_out = 0; 

Index = zeros ( 1 , dim_Ctr) ; 
nxtVec = nVec; 
for jt = l:dim_nVec 
if nVec(jt)==0 

jindex = jIndex+1; 
Index ( j Index) = jt; 
end % If nVec 
end % For jt 

% Index [checked and works] 
for jt = l:2"dim_Ctr 
for jtt = l:dim_Ctr 

nxtVec ( Index ( jtt ) ) = ctab ( jt, jtt ) ; 
end % For jtt 
% nxtVec [checked and works] 
y_out = y_out + Ffn (CaseCode, nxtVec) ; 
% jt -report = jt % * drop in production 

end % For jt 
end % of big block 

y_out = (-1 ) " sum (nVec) . * y_out; 
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(8) Ffn.m 



function [ z_out ] = Ff n (CaseX, nVecX) 

% \mathcal{F} - function for various cases. Casecodes [1 Cpt] [2 mu] [3 n] 

% [4 Cpt k m] [5 mu k m] [6 n k m] 

% Cpt is a real number or integer that points to an array Cvec of coefficients. 

% called by Tfn.m which supplies CaseX and nVecX. 

global alphaO dark vsq % supplied by muRun or by setUp .m 

alpha = ones (1,4); 

darkfac = 1; 

usq = 1-vsq; 

for kt = 1:4 

if nVecX(kt)==l 

alpha (kt) = alphaO(kt); 

darkfac = - darkfac .*( 1-dark (kt ));% might vectorize later 

end % If nVecX 
end % For kt 

w = alpha ( 1 ) . *alpha ( 4 ) . *usq; 

X = alpha ( 1 ). *vsq; 

y = alpha (2) . *alpha (3) . *usq; 

z = alpha (2 ) . *vsq; 

switch CaseX (1) 

case 1 

z_out = GCf n (w, X, Y, z ) ; % NOT YET IMPLEMENTED 
case 2 

mu_arg = CaseX (2); 
z_out = Gmuf n (mu_arg, w, X, y, z ) ; 
case 3 

n_arg = CaseX(2); 
z_out = Gnf n (n_arg, w, X, y, z ) ; 
case 4 

k_arg = CaseX(2); 

m_arg = CaseX(3); 
z_out = GCkmfn (k_arg,m_arg,w,x,y, z) ;% assumes global \mathbf{c} 
case 5 

mu_arg = CaseX (2); 

k_arg = CaseX (3); 

m_arg = CaseX (4); 
z_out = Gmu_kmf n (mu_arg, k_arg, m_arg, w, X, y , z ) ; 
case 6 

n_arg = CaseX (2) 

k_arg = CaseX (3) 

m_arg = CaseX (4) 

z_out = Gnkmf n (n_arg, k_arg, m_arg, w, X, y , z ] 
otherwise 

end 

z_out = z_out .* darkfac; 



91 



(9) Gkmfn.m 

f unction [y_out ] = Gnkmf n (n, k, m, w, x, y, z ) 
% 4 NOV 04 

% "vectorize" over w , x, y, z as well as over zetaX 
% works only if 32 \ge n \ge k + m 

% implemented to assume zeta = \infty if zeta \ge 10'~40 

% gamTab{n+l} = factorial (n) For \le n \le 32 (think gamma function). 

global XiAr CaseSetUp gamTab 
zetaX = CaseSetUp{5}; 
if zetaX==0 % BIG BLOCK 

aO = (gamTab{n-k+l}/ (gamTab{n-k-m+l}*gamTab{m+l}* (n+1) ) ) . . . 
.*x.''k.*z.''m.*y.'' (n-k-m) ; 
if n-k-m==0 

y_out = aO; 
else 

avec = ones ( 1 , n-k-m+1 ) ; 
avec ( 1 ) = aO ; 

for jt =1 : n-k-m 

avec(jt+l) = ( jt+k) * (n-k-m+l-jt) / ( (n-k+1- jt) * jt) ; 

end % For jt 
xt = w. /y; 

y_out = sumf n (avec, xt ) ; 
end % If n-k-m==0 
else if zetaX < 1 ^ 4 % General Case 
zlen = length ( zetaX) ; 
NsqFac = zeros ( zlen, 1 ) ; 
for jt = 0:n 

NsqFac = NsqFac + XiAr{ jt + l} . *XiAr{n- 
jt + 1} . / (gamTab{ jt + l}*gamTab{n- jt + l}) ; 
end % For jt 

NsqFac = (gamTab{n+l}) . *NsqFac; 
NsqFac = 1. /NsqFac; 

Const = NsqFac .* (gamTab{n+l}/ (gamTab{m+l}*gamTab{k+l}) )*.. . 
(x.'~k.*z.'~m.*y." (n-k-m) ) ; 

Tot = zeros ( zlen, length (w) ) ; 

for j = : n-k-m 

Tot = Tot + XiAr{ j+k+1} . *XiAr{n-k- j+1}* (w. /y) . " j . / . . . 
(gamTab{ j+l}*gamTab{n-k-m- j+1}) ; 
end % For j 
y_out = Tot.*Const; 
else % limit as zeta --> \infty 

y_out = 2" (-n) *gamTab{n+l} . *x. "k. *z . "m . * (w+y ) . " (n-k- 
m) ./. . . 

(gamTab{k+l}*gamTab{m+l}*gamTab{n-k-m+l}) ; 

end 

end % BIG BLOCK 
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(10) sumf n .m 

function [y] = sumfn(ar,x) 

% ar is a vector of the fomi [a_0, a_l/a_0, a_2/a_l, . . . , a_n/a_{n-l}] 
% X is a variable value or a vector of variable values 
%y = sum_{j = 0}"n a_j.*x.''j. 
dimar = size (ar) ; 

ntemp = dimar(2); %biggestnis ntemp - 1; 
% setup ztemp 
dimx = size (x) ; 

nxtemp = dimx (2); 
ztemp = ones ( 1 , nxtemp) ; 
for ct = ntemp : -1 : 2 % trouble if ntemp < 2. 
y = X . *ar (ct ) . *ztemp; 
ztemp = 1+y; 
end 

y = ar (1) . *ztemp; 

%testbyar = [2 1.2 1.2 1.2 1.2], x = [2 3 4] 
%sumfn([2 1.2 1.2 1.2 1.2], [2 3 4]) Checks. 
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(11) Gnfn.m 

function [y_Gn] = Gnf n (n_arg, w, x, y, z ) 

% 26 OCT 04 % does limit of z eta — > \infty if zetaX > 10''40 

% 29 OCT tested general zet aX against limits OK. 

% Old Gnf n ( zetaX, n_arg, w, x, y , z ) 

global XiAr CaseSetUp gamTab 

xl = (y+z) ; 

x2 = (w+x) ; 

zetaX = CaseSetUp{5}; 

if zetaX == 

if xl == x2 

y_Gn = xl."n_arg; 

else 

y_Gn = (xl . " (n_arg+l) -x2 . " (n_arg+l) ) . / ( (n_arg+l) . * (xl-x2) ) ; 

end % If xl 
else 

if zetaX < 10^40 

zlen = length ( zetaX) ; 

NsqFac = zeros ( zlen, 1 ) ; 
for jt = 0:n_arg 

NsqFac = NsqFac + XiAr{ jt + l} . *XiAr{n_arg- jt + l} . / . . . 

(gamTab{ jt + l}*gamTab{n_arg- jt + 1}) ; 
end % For jt 

NsqFac = (gamTab{n_arg+l}) . *NsqFac; 

NsqFac = 1. /NsqFac; 

Const = NsqFac* (y+z ). '~n_arg; 

% Modified from Gnkmf n 

Tot = zeros ( zlen, length (w) ) ; 
for j = : n_arg 

Tot = Tot + XiAr{ j+1} . *XiAr{n_arg- j+1}* . . . 
( (w+x) . / (y+z) ) . " j . / (gamTab{j+l}*gamTab{n_arg- j+1}) ; 

%(*) 

end % For j 

Tot = gamTab{n_arg+l} . *Tot ; 
y_Gn = Tot.*Const; 
else % limit as zeta --> \infty 

y_Gn = ( (xl+x2 ) . /2 ) . "n_arg; 
end % If zetaX < 10^40 
end 
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(12) Gmu_kmfn.m 



function [y_ans ] = Gmu_kmf n (mu_arg, k_arg, m_arg, w, x, y, z ) 
% 30 OCT 04 

global CaseSetUp gamTab 

zetaX = CaseSetUp{5}; 

% checked sum against Gmuf n, OK. 

if zetaX < 10^40 %zeta \ge 1 M treated as infinite 
nmax = floor (3*mu_arg+16) ; 
y_ans = 0; 

for nt = k_arg + m_arg:nmax 
incr = Gnkmf n (nt , k_arg, m_arg, w, X, y , z ) ; 

y_ans = y_ans + mu_arg"nt . *incr . /gamTab{nt+l}; 
end % For nt 
y_ans = exp(-mu_arg) .*y_ans; 
else % limit as zeta — > \infty 
mu2 = mu_arg/2; 

y_ans = exp (-mu2 . * (2-w-y) ) . * (mu2 . *x) . ''k_arg. * (mu2 . *z) . "m_arg. / . . . 
(gamTab{k_arg+l}*gamTab{m_arg+l}) ; 

end 
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(13) Gmuf n .m 

function [y_out ] = Gmuf n (mu_arg, w, x, y, z ) ; 

global CaseSetUp gamTab 

% gamTab{n+l} = f actorial{n} For n = 0, nmax 
zetaX = CaseSetUp{5}; 
if zetaX == 
xl = (y+z) ; 
x2 = (w+x) ; 
if mu_arg == 

y_out = 1; 
else 

if xl == x2 

y_out = exp (-mu_arg . * ( 1-xl ) ) ; 
else 

y_out = exp (-mu_arg) . * (exp (mu_arg . *xl ) 

-exp (mu_arg . *x2 ))./... (mu_arg . * (xl-x2 ) ) ; 
end % If xl == x2 
end % If mu_arg == 
else if zetaX > 10^40 

y_out = exp (-mu_arg . * (2-w-x-y-z ) . /2 ) ; 
else % General Case < zetaX \le 10^40 

y_out = zeros ( length ( zetaX) , length (w) ) ; 
% column vector of same length as z e t aX 
nmax = f loor (3*mu_arg+16) ; 
for nt = 0:nmax 

incr = Gnf n (nt , w, X, y , z ) ; 

y_out = y_out + mu_arg"nt . *incr . /gamTab{nt+l}; 
% gamTab{nt+l} = factorial (nt ) 
end % For nt 

y_out = exp (-mu_arg) . *y_out ; 
end % General Case 
end 
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